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Abstract 

We  generalize  the  hydrodynamic  lattice  gas  model  to  include  arbitrary  numbers  of  particles  moving 
in  each  lattice  direction.  For  this  generalization  we  derive  the  equilibrium  distribution  function  and 
the  hydrodynamic  equations,  including  the  equation  of  state  and  the  prefactor  of  the  inertial  term  that 
arises  from  the  breaking  of  galilean  invariance  in  these  models.  We  show  that  this  prefactor  can  be  set  to 
unity  in  the  generalized  model,  therby  effectively  restoring  galilean  invariance.  Moreover,  we  derive  an 
expression  for  the  kinematic  viscosity,  and  show  that  it  tends  to  decrease  with  the  maximum  number  of 
particles  allowed  in  each  direction,  so  that  higher  Reynolds  numbers  may  be  achieved.  Finally,  we  derive 
expressions  for  the  statistical  noise  and  the  Boltzmann  entropy  of  these  models. 


20050829  018 


1 


»• 


1  Lattice  Gases 

Lattice  gas  automata  (LGA)  are  a  class  of  dynamical  systems  in  which  particles  move  on  a  lattice  in  discrete 
time  steps.  If  the  collisions  between  the  particles  conserve  mass  and  momentum,  and  the  lattice  is  sufficiently 
symmetric,  the  coarse-grained  behavior  of  the  system  can  be  shown  to  be  that  of  a  viscous  fluid  in  the 
appropriate  scaling  limit  [1,  2,  3,  4],  Used  as  an  algorithm  for  simulating  hydrodynamics,  the  method  has 
the  virtues  of  exact  conservation  laws,  and  of  unconditional  numerical  stability. 

In  a  typical  LGA,  there  is  an  association  between  the  lattice  vectors  and  the  particles  at  each  site.  If 
there  are  n  lattice  vectors,  then  the  state  of  the  site  is  represented  by  n  bits.  Each  bit  represents  the  presence 
or  absence  of  a  particle  in  the  corresponding  direction.  At  each  time  step,  a  particle  propagates  along  its 
corresponding  lattice  vector  and  then  collides  with  other  arriving  particles  at  the  new  site.  (Note  that  rest 
particles  can  be  subsumed  into  this  scheme  by  associating  them  with  null  lattice  vectors.)  The  collisions  are 
required  to  conserve  particle  mass  and  momentum. 

Relevant  dimensionless  quantities  of  a  LGA  are  the  Knudsen  number,  Kn,  defined  as  the  ratio  of  the 
mean-free  path  to  the  characteristic  length  scale;  the  Strouhal  number,  Sh,  defined  as  the  ratio  of  the  mean- 
free  time  to  the  characteristic  time  scale;  the  Mach  number,  M,  defined  as  the  ratio  of  the  characteristic 
velocity  to  the  speed  of  sound;  the  Reynolds  number.  Re  ~  M/Kn;  and  the  fractional  variation  of  density 
from  its  average  value,  5p/p.  Hydrodynamic  behavior  [5]  is  attained  in  the  limit  as  Kn  and  Sh  go  to  zero. 
Viscous  hydrodynamics  [5]  is  attained  when  Sh  ~  Kn^  in  this  limit.  Incompressible  viscous  hydrodynamics  [6] 
is  then  attained  when  we  also  have  M  ~  Kn  so  that  Re  ~  0(1),  and  Sp/p  ~  Kn^. 

The  Chapman-Enskog  procedure  is  a  perturbation  expansion  in  the  above-described  asymptotic  ordering. 
For  a  LGA  whose  collisions  conserve  mass  and  momentum  on  a  lattice  of  sufficient  symmetry  (quantified 
below),  the  local  equilibrium  distribution  function  can  be  shown  to  be  Fermi-Dirac  in  nature  [2,  3,  4].  The 
Chapman-Enskog  procedure  can  then  be  used  to  compute  the  correction  to  this  Fermi-Dirac  distribution  and 
thereby  show  [1,  2,  3,  4]  that  the  pressure,  P,  and  the  momentum  density,  u,  obey  the  following  equations 
in  the  asymptotic  limit: 

Vu  =  0 

5  +  — u  ■  Vu  =  -VP  +  u{p)V^u, 

Ot  p 

where  p  is  the  fluid  density  (a  constant  in  this  limit).  The  analysis  also  yields  expressions  for  the  functions 
g{p)  and  v{p),  and  an  equation  of  state  for  P.  In  particular,  the  form  of  these  equations,  the  equation 
of  state,  and  the  expression  for  the  function  g{p)  depend  only  on  the  fact  that  mass  and  momentum  are 
conserved  -  and  are  the  only  things  conserved  -  by  the  collisions.  The  expression  for  the  viscosity,  u{p) 
depends  on  the  details  of  the  collision  rules  used. 

Since  the  fluid  density  p  is  a  constant  in  the  asymptotic  limit,  the  factors  g{p)  and  v{p)  are  also  constants. 
As  has  been  noted,  the  latter  is  the  fluid  viscosity.  The  presence  of  the  former  is  reflective  of  a  breaking  of 
galilean  invariance,  due  to  the  fact  that  the  lattice  itself  constitutes  a  preferred  galilean  frame  of  reference. 
For  a  single-phase  LGA,  the  former  factor  can  eetsily  be  scaled  away  by  redefining  the  momentum  density 
and  pressure  as 

U  =  g{p)vi 

and 

V  =  g{p)P, 

where  u  and  P  are  those  measured  in  the  simulation.  For  compressible  flow,  or  for  multiphase  flow  with 
interfaces,  however,  the  presence  of  the  g(p)  factor  is  problematic,  and  various  techniques  have  been  proposed 
to  remove  it.  It  has  been  shown  that  this  can  be  done  by  judiciously  violating  semi-detailed  balance  in  the 
collision  rule  [7],  or  by  adding  many  rest  particles  at  each  site  [8]. 

The  unconditional  stability  of  the  lattice  gas  procedure  arises  from  a  requirement  that  the  collisions 
satisfy  a  statistical  reversibility  condition  known  as  semi-detailed  balance  (SDB).  The  collision  process  is 
fully  specified  by  the  transition  matrix  A(s  — >  s')  which  is  the  probability  that  the  incoming  state  s  will 


2 


result  in  outgoing  state  s'.  Since  collisions  must  result  in  some  outgoing  state,  conservation  of  probability 
requires  that 

=  (1.1) 

s' 

SDB  is  then  the  condition  that 

^yl(s^s')  =  l.  (1.2) 

a 

(Note  that  the  condition  of  detailed  balance  (DB),  j4(s  — >  s')  =  i4(s'  — )•  s),  implies  that  of  SDB,  but  not  vice 
versa;  that  is,  SDB  is  a  weaker  condition  than  DB.)  FVom  SDB,  it  is  possible  to  prove  an  /f-theorem,  from 
which  follows  the  unconditional  stability  of  the  lattice  gas  algorithm. 

An  important  limitation  of  the  lattice  gas  procedure  has  to  do  with  the  statistical  noise  associated  with 
the  coarse-grained  averaging  that  is  necessary  to  get  the  hydrodynamic  quantities  that  obey  the  above 
fluid  equations.  For  n  bits  per  site,  and  for  coarse-grained  averages  over  blocks  of  N  sites,  the  noise  is 
of  order  ~  Ijs/nN.  For  some  applications  -  most  notably  the  simulation  of  complex  fluids  -  a  certain 
controllable  amount  of  noise  is  actually  desirable  because  it  is  essential  to  the  physics;  for  simple  fluid 
dynamics  computations,  on  the  other  hand,  the  noise  is  a  nuisance. 


2  Lattice  Boltzmann  Equations 

Because  of  their  noise  and  lack  of  galilean  invariance,  LGA  have  been  replaced  by  Lattice  Boltzmann  Equa¬ 
tions  (LBE)  for  many  hydrodynamics  applications  of  interest  in  recent  years  [9] .  These  methods  keep  track 
only  of  an  averaged  occupation  number  of  particles  in  each  direction  at  each  site.  Moreover,  the  collision 
operator  most  often  used  is  a  simple  relaxation  to  a  noiseless  equilibrium,  thereby  eliminating  the  statistical 
fluctuations  that  are  inherent  in  the  LGA  method.  This  means  that  in  complex  fluid  applications  for  which 
statistical  fluctuations  are  an  essential  part  of  the  physics,  they  have  to  be  reintroduced  artificially  [10]. 

For  a  lattice  Boltzmann  equation  corresponding  to  a  lattice  gas  with  only  one  bit  per  lattice  vector,  this 
real- valued  distribution  function  is  bounded  between  zero  and  one.  This  need  not  be  the  case,  however,  and 
the  LBE  procedure  allows  one  to  tailor  the  equilibrium  distribution  function  to  satisfy  certain  desiderata. 
Among  these  is  the  ability  to  demand  galilean  invariance  (g{p)  =  1)  [11]. 

At  the  same  time,  the  LBE  method  gives  up  two  of  the  principal  advantages  of  LGA’s:  Due  to  the  roundoff 
error  inherent  in  manipulations  of  real  numbers  on  a  computer,  it  no  longer  maintains  the  conservation  laws 
exactly.  Moreover,  LBE’s  are  no  longer  unconditionally  stable;  indeed,  they  are  subject  to  a  variety  of 
numerical  instabilities,  most  of  which  are  not  well  understood. 

3  Integer  Lattice  Gas  Automata 

In  this  paper,  we  investigate  a  simple  generalization  of  the  lattice  gas  concept  that  can  be  used  to  control 
the  level  of  statistical  fluctuations  -  reducing  it  if  desired,  but  not  necessarily  eliminating  it  altogether  - 
while  maintaining  the  conservation  laws  exactly,  preserving  unconditional  stability,  and  allowing  for  galilean 
invariance. 

The  use  of  a  single  bit  per  each  of  n  directions  to  represent  the  state  of  a  given  lattice  site  means  that 
the  number  of  particles  moving  along  any  lattice  direction  is  either  zero  or  one.  We  generalize  this  by 
allowing  for  up  to  L  bits  per  direction,  for  a  total  of  nL  bits  per  site,  so  that  the  number  of  particles  moving 
along  any  lattice  direction  can  range  from  0  to  2^  —  1.  The  total  number  of  states  per  site  is  then  2”^. 
Computationally,  this  means  that  the  state  of  each  direction  is  described  by  an  integer  of  L  bits;  hence,  the 
terminology.  Integer  Lattice  Gas  Automata  (ILGA). 

To  simplify  the  derivation  of  the  hydrodynamic  equations  of  an  ILGA,  we  use  the  Boltzmann  molecular 
chaos  approximation,  so  that  all  quantities  in  our  analysis  are  ensemble  averaged,  and  we  indiscriminately 
commute  the  application  of  this  average  with  the  collision  process.  We  also  assume  that  the  particles  are  of 
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unit  mass.  Denote  the  ensemble-averaged  value  of  the  ^th  bit  in  the  I'th  direction  by  where  0  <  i  <  n—  1 
and  0  <  ^  <  L  —  1.  Also,  denote  the  lattice  vector  for  the  zth  direction  by  Cj.  The  distribution  function  for 
the  total  number  of  particles  in  each  direction  is  then, 

L-l 

N'  =  Yl  (3.3) 

e=o 


The  ensemble-averaged  mass  and  momentum  densities  are  then  given  by 

n— 1  n—lL—l 


(3.4) 


and 


n— 1 


t=0  ^=0 


n-1 L-l 


=  E  E  (3-5) 

t=o  >=o  r=o 

Let  us  also  associate  an  energy  Si  with  each  particle  in  direction  i.  The  ensemble-averaged  energy  density  is 
then  given  by 


n— 1 


n-1 L-l 


(3.6) 


t=0 


i=0  ^=0 


4  Thermodynamics  of  the  Integer  Lattice  Gas 

We  first  consider  the  thermodynamics  of  the  integer  lattice  gas.  The  grand  canonical  partition  function  is 

Z=Y.^xp[-p{E-a■P-^lM)], 

{N] 

where  the  sum  is  over  all  possible  states  of  the  ILGA  (that  is,  each  A'(x)  is  summed  from  0  to  2^  —  1), 
where  /9,  a  and  fj,  are  Lagrange  multipliers,  and  where 

V  n 

m^eE^'w- 


X  i 

are  the  total  mass,  momentum  and  energy,  respectively,  of  all  the  particles  on  a  lattice  of  V  sites.  Thus,  we 
have 


Z  = 


Eexp 

{JV} 


V  n 

ij.)  A‘(x) 

X  i 


E  nn®^P  [-/3(e<  -a-Ci-fj.)  iV'(x)] 

{AT}  X  i 
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V  n  2''-! 


V  n  2''-! 


nn  exp  [-^ (Ci  -  a  •  Ci  - /i)  fc]  =  JJ  JJ 

X  i  fc=0  X  i  k—0 


where  we  have  defined  the  fugacity 


The  grand  potential  is  then 


■i-a 

n 

1 

n 

Ein 

i 

A  -{2^2- 

l-2< 

yi 

2^(2^^^  \ 

- 

1  -  (2^2^  ) 

J' 

so  that 

5^  (^fl)  =  -  5^  In  2  -  V  ^  j  (s,  ^  .  c,  -  rt . 

We  identify  the  equilibrium  distribution  function 

which  gives  the  mean  number  of  particles  moving  in  each  lattice  direction.  Since  this  has  a  maximum  of 
2^  —  1 ,  we  also  define  the  fractional  occupation  number 


hiz)  = 


_  FUz) 


2^-r 


Figure  1  shows  fdz)  plotted  against  2  for  several  values  of  L. 

In  terms  of  the  equilibrium  distribution  function,  we  have 

n  +  f3^  =  n-T^^vf^FUz')  {Bi-a-a-fz), 

^  i 

where  T  =  1//?  is  the  temperature.  It  follows  that 

n  =  {H)-a-{P)-,,{M)-T{S), 
where  we  have  identified  the  average  energy 

n 

{H)  =  Vj2FLiz')ei, 

i 

the  average  momentum 

{p)  =  vf;^FUzicu 

i 

the  average  mass 

(M}  =  V'£Fl{z% 


the  average  mass 
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Fractional  Occupation  Number 


Figure  1;  fi{z)  versus  z  for  several  values  of  L.  The  black  curves  represent  L  values  from  1  to  6,  with 
increasing  steepness,  while  the  gray  curve  is  the  limit  as  L  — >  oo. 
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Figure  2;  Entropy  excess  ASi,  versus  fractional  occupation  number  fL{z).  The  black  curves  represent  L 
values  from  1  to  6,  increasing  downward,  while  the  gray  curve  is  the  limit  as  L  oo. 

and  the  average  entropy 

i 

In  the  expression  for  the  entropy  we  have  defined  the  function 

Sl{z)  =  In  (l  -  j  In  -  In  (1  -  z)  -  In  W  (4-9) 

as  the  entropy  per  lattice  direction.  Thus,  in  addition  to  the  form  for  the  equilibrium  distribution  function, 
this  analysis  has  provided  us  with  an  expression  for  the  entropy  that  is  additive  in  the  contributions  from  each 
lattice  direction.  In  fact,  it  is  straightforward  to  show  that  — >•  L  In  2  in  the  limit  of  large  L,  corresponding 
to  a  dominant  contribution  of  In  2  per  bit  of  state.  The  excess 

A5i  =  5i,-Lln2  (4.10) 

is  then  0{1)  in  L  and  is  plotted  against  the  fractional  occupation  number  fi{z)  in  Fig.  2.  This  can  be 
interpreted  as  indicating  that  the  bits  are  most  random  at  half  filling;  elsewhere,  the  entropy  is  lower  than 
L  In  2  per  bit. 
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5  Kinetic-Theoretical  Treatment 


As  an  alternative  to  the  preceding  thermodynamic  treatment  of  the  integer  lattice  gas,  we  can  derive  the 
principal  results  from  a  kinetic-theoretical  argument.  For  example,  to  derive  the  form  of  the  equilibrium 
distribution  function,  Eq.  (4.8),  we  can  note  that  the  equilibrium  distribution  function  for  each  bit  must  still 
be  Fermi-Dirac  in  form,  since  each  individual  bit  is  either  occupied  or  not.  Thus, 


°  1  -b  exp  [2^/3  {si-a  ci-  /r)]  ’ 

where  the  multipliers  /3,  a  and  fj,  are  determined  in  terms  of  the  mass,  momentum  and  energy  densities  by 
their  definitions,  Eqs.  (3.4),  (3.5)  and  (3.6).  In  terms  of  the  fugacity,  Eq.  (4.7),  the  above  may  be  written. 


(5.11) 


The  equilibrium  distribution  function  for  each  direction  is  then  given  by  Eqs.  (3.3)  and  (5.11), 


L-i 


t=o 


where  we  have  defined  the  function 


We  see  that  this  series  telescopes  to  yield  the  form  derived  in  the  previous  section, 

nL 


(5.12) 


(5.13) 


6  Form  of  the  Hydrodynamic  Equations 

To  derive  the  hydrodynamic  equations,  we  first  expand  the  equilibrium  distribution  function  in  the  Mach 
number.  Here  and  henceforth,  we  specialize  to  the  case  of  no  internal  energy,  so  that  £;  =  0,  and  we  can 
absorb  the  multiplier  (3  into  a  and  /r.  TVeating  a  as  a  small  quantity,  the  fugacity  can  be  written 

=  e**  ^1  -b  a  •  Ci  ^aa  :  =  Zq  +  z\  +  z^, 

where  the  subscripts  of 

20  =  e'^ 

z\  =  zoa  •  a 

i  _  ^0 

22  =  ;  CjCj 

denote  the  order  of  the  Mach  number  expansion,  and  we  note  that  zq  is  independent  of  the  direction  i.  It 
follows  that 

N^  =  FUz')  =  FL{zo  +  z{+zi). 

Taylor  expanding,  we  get 

A'o  =  •^£,(•^0)  +  zoFi^{zo)a  ■  c*  -b  ^zo  N-FK^o)]^  aa  :  CjCi. 
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To  proceed,  we  must  make  some  assumptions  about  the  symmetries  of  the  lattice.  We  demand  that 


(6.14) 

i=0 

for  0  <  fc  <  4,  where  ®  denotes  the  outer  product,  and  where  is  the  completely  symmetric  and  isotropic 
tensor  of  rank  k, 

lo=l 
(li)i  =  0 
{l2)y  =  Sij 
(l3)yfc  =  0 

—  ^ij^kl  T  ^ik^jl  “t  ^il^jk‘ 

Note  that  Eq.  (6.14)  defines  the  coefficients  Ak  for  a  given  lattice. 

We  now  demand  that  ^ 

P=^Ni  =  AeFaze)  +  ^zo[zoF'M^;, 
i=o 


T*  — A 

u  =  ^  CiNo  =  A2ZoFlizo)a. 


If  we  now  let  z  denote  the  solution  to  the  equation 

^=Fx(^),  (6.15) 

it  follows  that  the  difference  between  z  and  zq  is  of  second  order  in  the  Mach  number,  so  that  we  can  solve 
for  u,  and  a.  We  find  that 

u 

A^zFiizy 

and  that  /x  =  In  zo  where  zq  is  the  solution  to  the  equation 

f 

Ao  2AoA2[zFl{z)f 

Inserting  these  results  into  the  distribution  function,  we  find 


“  >lo  A2  2AUzFiiz)] 

where,  again,  z  is  defined  by  Fi{z)  =  p/Ao- 

The  inviscid  part  of  the  stress  tensor  is  then  given  by 


l[^F’Az)f\'^  AoV’ 


A  Az  zlzFAz)]'  ,  A 

y  CiCiNn  =  -e/9  +  — - —  9  -7^14  -  -r^l2®  I2  : 

^0  2A2mz)f\A2^  Ao  "  V 

^  'A2  z[zFl{z)]'  (A,  A2\  2I 

Ao  ^  2 A2  [zF'j^  (z:) jH  ^2  Ao  /  ^ 


i.,^2 

=  -P(P>m)12  +  P(/2)  — , 
P 


+  AfzF'dz)r^^ 
Ai  [zFi{z)y 
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where  we  have  identified  the  factor  that  multiplies  the  inertial  term  in  the  Navier-Stokes  equations, 


and  the  equation  of  state, 


_  AoA,zFaz)mz)]' 

Al[zFi{z)f  ’ 


Eqs.  (6.15),  (6.17)  and  (6.18)  are  the  principal  results  of  this  section.  Eq.  (6.15)  gives  p  in  terms  of  the 
parameter  z.  Eq.  (6.17)  then  gives  g  in  terms  of  z,  so  that  Eqs.  (6.15)  and  (6.17)  are  a  pair  of  parametric 
algebraic  equations  for  g  in  terms  of  the  density  p.  Finally,  Eq.  (6.18)  gives  the  equation  of  state  for  P  in 
terms  of  p  and  u.  The  coefficients  Aj  that  appear  in  these  equations  are  given  in  terms  of  the  lattice  vectors 
by  the  conditions,  Eq.  (6.14). 

7  Example:  Bravais  Lattice 

As  a  concrete  example  of  this  formalism,  we  consider  the  case  of  a  regular  Bravais  lattice.  Examples  of  such 
lattices  with  the  requisite  symmetry  conditions,  Elq.  (6.14),  are  the  triangular  lattice  in  two  dimensions  [1] 
and  the  face-centered  hypercubic  lattice  in  four  dimensions  [3].  In  addition  to  the  n  directions  corresponding 
to  unit-speed  particles,  we  include  Ur  null  lattice  vectors  to  accomodate  rest  particles.  In  this  situation, 

Ao  —  n  +  Ur 


^  D(D  +  2)  ’ 

where  D  is  the  number  of  dimensions.  Inserting  these  into  Eqs.  (6.15)  through  (6.18),  we  find 

p  =  {n  +  nr)FL{z), 


Here  we  have  defined  the  function. 


which  we  plot  against  the  fractional  occupation  number. 


hiz)  = 


(2^  —  1)  (n  -1-  Ur)  2^  —  1  ’ 


for  several  different  values  of  L  in  Fig.  3.  For  L  =  1  it  is  a  straightforward  exercise  to  show  that  we  recover 
the  well  known  result  [3], 

which  decreases  monotonically  from  unity  at  /  =  0,  to  zero  at  half-filling  (/  =  1/2),  after  which  it  becomes 
negative.  For  L  >  1,  we  see  that  this  decrease  is  no  longer  monotonic,  since  the  slope  at  the  origin,  ^'(0), 
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is  positive.  Thus,  for  L  >  1,  the  function  g  has  a  maximum  for  some  0  <  /  <  1/2.  The  location  of 
this  maximum  approaches  /  =  0  as  Z,  — >  oo.  (The  limit  of  infinite  integers,  i.e.,  L  — >■  oo  is  discussed  in 
Appendix  A,  and  is  shown  as  a  shaded  curve  in  Fig.  3.) 

Galilean  invariance  is  achieved  when  p  =  1,  or 


If  the  quantity  (1  +  2ID)nl{n  +  n^)  is  greater  than  the  maximum  value  of  G^,  then  galilean  invariance  is 
impossible  for  those  values  of  D,  n,  rir,  and  L;  if  it  is  less  than  this  maximum,  then  there  are  two  densities 
at  which  galilean  invariance  is  achieved.  Some  of  these  values  are  tabulated  for  the  FHP  and  FCHC  lattice 
gases  in  Fig.  4. 


8  Viscosity 

To  compute  the  viscosity  of  a  ILGA  in  the  Boltzmann  molecular  chaos  approximation  [3],  we  consider  its 
ensemble-averaged  collision  operator,  This  quantity  is  the  ensemble  average  of  the  increase  in  bit  £  in 
direction  i  due  to  collisions.  It  is  given  by 

a, a' 


where  'P{s)  is  the  probability  that  the  incoming  state  is  s,  A(s  — >■  s')  is  the  probability  that  the  collision 
process  takes  incoming  state  s  to  outgoing  state  $',  and  is  the  value  of  bit  £  in  direction  i  in  incoming 
state  s  (and  likewise  for  outgoing  state  s').  In  the  Boltzmann  approximation,  the  probability  of  a  state  s  is 
the  product  of  the  corresponding  fractional  occupation  numbers,  or  their  complements. 


v=oy=o 

To  get  the  total  increase  of  particles  in  direction  t,  we  take  the  sum 

£.-1  _  n— 1  L-\  k' 

fl*  =  ^  ^  A(s  s')(s'*  —  .<)  n  n  ("‘'0'  (■-«*''')■■ 

^=0  s,s*  Jfc'=0/=0 


s‘  =  ^  2V-' 


is  the  total  number  of  particles  in  direction  i  in  state  s  (and  likewise  for  s'). 

To  compute  the  viscosity,  we  must  form  the  Jacobian  matrix  of  the  collision  operator.  Direct  calculation 
yields 

dJSfk,]  'iv*.j(l  -  A/’fc.j)’ 

We  would  like  to  evaluate  this  Jacobian  at  the  equilibrium  given  by  Eq.  (5.11), 

^ 

°  1  +  (^*)-2^  ■ 

Taking  the  derivative  of  this  equation  with  respect  to  the  fugacity. 


2^(z^)-2-i 

[1  +  {zk)-vf  zfc 
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FHP  Lattice  Gas  (Z)  =  2,  n  =  6) 

IB 

Low-Density  Root  High-Density  Root 

o 

0.0 

0.0 

1 

6 

0.0396831 

0.143848 

oo 

0.0 

0.168451 

2 

4 

0.0704358 

0.126560 

5 

0.0322583 

0.177356 

6 

0.0158730 

0.195949 

CXD 

0.0 

0.212636 

4 

3 

0.0362392 

0.167555 

4 

0.0166667 

0.228582 

5 

0.0080645 

0.253770 

6 

0.0039683 

0.265426 

OO 

0.0 

0.276535 

FCHC  Lattice  Gas  (D 

=  4,  n  =  24) 

IB 

IB 

Low-Density  Root 

High-Density  Root 

0 

4 

0.0704358 

0.126560 

5 

0.0322583 

0.177356 

6 

0.0158730 

0.195949 

OO 

0.0 

0.212636 

1 

4 

0.0528917 

0.152419 

5 

0.0253456 

0.193030 

6 

0.0124717 

0.209795 

OO 

0.0 

0.225163 

2 

4 

0.0417410 

0.171769 

5 

0.0201613 

0.207212 

6 

0.0099206 

0.222576 

OO 

0.0 

0.236849 

4 

3 

0.0654937 

0.120009 

4 

0.0266677 

0.203001 

5 

0.0129032 

0.232239 

6 

0.0063492 

0.245488 

OO 

0.0 

0.257994 

Figure  4:  Values  of  /  G  (0, 1/2)  such  that  g  =  1 
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we  can  use  the  chain  rule  to  get  the  integer  version  of  the  Jacobian  of  the  collision  operator  at  equilibrium, 

an' 

0 

an'-^  I  aN'^'^az^  \ 


aN>‘ 

L-l 

e,]=o 


0  am/az'' 


=  E-*(«  ^  . 


(8.20) 


In  fact,  we  need  this  result  only  in  the  limit  of  zero  Mach  number,  so  we  can  use  the  lowest  order  expression 
for  the  fugacity,  z*'  =  z  (see  Sec.  6),  which  is  independent  of  the  index  k.  We  find  that  the  zero  Mach  number 
limit  of  the  Boltzmann  probability  of  state  s  is  given  by 


fc'=o j'=o  y=o  /  \i-r2 


-Py  («) 


where 


n— 1 


*=0 

is  the  total  number  of  populated  bits  in  the  jth  binary  digit.  It  follows  that 

rt-i 


PW- 


n('+^=T‘ 

J=0 


3=0 


,pU) 


where 


L-l 

P(«)  -  E 

3=0 

is  the  total  number  of  particles  present  in  state  s. 

Inserting  this  result  into  the  expression,  Eq.  (8.20),  for  the  collision  operator,  we  obtain 

=  7^)  (t^)"  E"^^"  ^ 

As  a  consequence  of  conservation  of  probability,  Eq.  (1.1),  and  semidetailed  balance,  Eq.  (1.2),  it  follows 
that  the  second  term  in  square  brackets  vanishes,  so  we  finally  get 


At  first  order  in  Knudsen  number,  the  kinetic  equation  is  [4] 

Ci  •  VAT*  =  J'jN(, 


where  there  is  an  understood  summation  over  j.  The  only  part  of  the  left-hand  side  that  contributes  to  the 
viscosity  comes  from  the  second  term  on  the  right-hand  side  of  Eq.  (6.16),  whence 


Now,  J  is  a  singular  matrix;  it  has  a  null  eigenvector  corresponding  to  each  hydrodynamic  mode  of  the 
system.  These  null  eigenvectors  span  what  we  shall  call  the  hydrodynamic  subspace  of  the  system.  The 
complement  of  this  subspace  is  called  the  kinetic  subspace,  and  is  spanned  by  the  kinetic  modes  with  nonzero 
(negative)  eigenvalue.  If  we  restrict  our  attention  to  the  kinetic  subspace,  then  we  can  form  the  pseudoinverse 
of  J,  denoted  by  J“^,  in  terms  of  which  we  may  write 

:  Vu. 

The  conservation  law  for  momentum  then  contains  the  term 


^  {aa  ■  VAT-  +  iciCiCi  :  VVJV*  +  •  •  •) 

t 

( 

1 


=  V- 


A2 


Y^CiCi{j  ^CjCj  +  l^CiCiCiCi 


:  (Vu)  \  +  ■ 


We  note  that  J  ^  is  diagonalized  and  degenerate  in  the  subspace  spanned  by  the  n  outer  products  of  the 
lattice  vectors  with  themselves;  that  is 


(8.21) 


where  A  is  a  scalar,  whence  the  above  term  in  the  momentum  conservation  equation  becomes 

[-44 


V- 


A2 


(-A+i)l4:(Vu) 


+  ...  =  V. 


A2 


(-A+l)Vu 


from  which  we  identify  the  kinematic  viscosity. 


The  quantity  A  is  then  determined  by  taking  the  double  spatial  dot  product  of  CjCj  on  both  sides  of  Eq.  (8.21), 
and  summing  over  i  to  get 

n  =  -A  ^  J'j  {a  ■  Cjf , 

whence 

s')  {s''  -  s')  id  .  . 

a^s' 

where  /  =  fii^)  determines  the  parameter  z  in  terms  of  the  fractional  occupation  number.  This  result  is 
easily  seen  to  reduce  to  that  of  Henon  [12]  when  L  =  1. 

We  computed  the  viscosity  of  an  L  =  2  lattice  gas  in  two  dimensions  {D  =  2)  by  measuring  the  decay 
of  a  shear  wave  in  periodic  geometry.  We  used  a  lattice  of  size  512  x  512  on  a  CAM-8  Cellular- Automata 
Machine  [13].  The  probabilistic  collision  procedure  used  obeyed  semi-detailed  balance,  with  each  outgoing 
state  allowed  by  the  conservation  laws  sampled  with  equal  probability.  Fig.  5  shows  the  decay  of  the  shear 
wave  amplitude  to  be  exponential  in  nature,  as  is  appropriate  for  Navier-Stokes  evolution.  The  time  constant 
of  the  exponential  then  determines  the  viscosity,  which  is  plotted  as  a  function  of  density  in  Fig.  6,  along  with 
the  curve  predicted  by  the  theory  given  above.  While  the  agreement  with  theory  is  good  at  intermediate 
values  of  the  fractional  occupation  number  near  half  filling,  we  note  that  it  is  seriously  in  error  at  low 
(and  high)  fractional  occupation  numbers.  At  present,  we  attribute  this  discrepency  to  deviations  from  the 
Boltzmann  molecular  chaos  approximation,  and  we  plan  to  investigate  them  using  kinetic  ring  theory  [4]  in 
a  forthcoming  publication. 


-1 

A  nzF'^{z)[l-z^‘-J 
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Kinematic  Viscosity 


9  Statistical  Noise 


Finally,  we  consider  the  statistical  noise  of  the  ILGA  model.  With  the  maximum  number  of  particles  per 
direction  increasing  as  2^,  one  might  naively  expect  the  noise  level  to  decrease  with  L  as  \l\p2p  ~ 
Unfortunately,  as  we  shall  show,  this  expectation  is  not  realized,  due  to  the  extremely  narrow  dynamic  range 
of  the  fugacity  for  large  L.  This  is  best  seen  in  Fig.  1,  in  which  the  effective  width  of  the  function  near 
z=\  decreases  like  2“^,  making  for  a  subtle  limiting  process  that  is  discussed  in  Appendix  A. 

Let  n*’*(x,t)  be  the  precise  value  of  bit  i  in  direction  i  at  lattice  site  x  at  time  t.  The  ensemble-average 
of  this  quantity  is  as  used  in  the  text  of  the  paper.  The  mean  number  of  particles  in  a  (space-time) 
block  of  N  sites  is  then  ^  ^ 

^^  =  EEE  2‘(n''\^,t))  =  nNFUz), 

{x,t)  i  ^=0 

where  the  angle  brackets  denote  the  ensemble  average. 

The  mean  square  of  the  number  of  particles  in  this  block  of  sites  is  then 

N  N  n  n  L~1  L-1 
(x,t)(x',t')  t  i'  e=oe'=o 


The  bits  are  either  zero  or  one,  and  in  the  Boltzmann  molecular  chaos  approximation  different  bits  are 
uncorrelated.  It  follows  that 


^n‘'^(x,t)n‘'’^'(x',f')^  =  (n‘’^(x,t))  ^n‘'’^'(x',t'))  +  (n*’^(x,t))  (l  -  (n‘’^(x,t)))  , 


whence 

=  JF2  +  nN  ^  ^  =  7?  -b  nNzFliz). 

It  follows  that  the  standard  deviation  of  the  number  of  particles  in  the  block  is  —  F\.  To  define  a 
fractional  noise,  we  could  divide  this  by  the  mean  number  of  particles,  but  it  preserves  particle- hole 
symmetry  if  we  instead  divide  it  by  the  square  root  of  the  product  of  the  mean  number  of  particles  and  the 
mean  number  of  holes,  thus 


F2-FI 


Fx\nN{2^-\)-Fx\ 


1  I  zf  (z) 

y/nN  (2^-1)]]  hiz)  [1  -  hiz)]  ■ 


(9.22) 


This  appears  to  decrease  exponentially  with  L,  but  it  must  be  noted  that  the  logarithmic  derivative  of  fi,{z) 
goes  as  2^  at  2  =  5.  Since,  for  fixed  fractional  occupation  number  fi,  z  tends  to  5  as  L  tends  to  infinity, 
we  see  that  AF  is  order  unity  in  L.  Thus,  the  fractional  noise  does  decrease  with  L,  but  not  as  rapidly  as 
one  might  hope.  It  is  plotted  as  a  fraction  of  l/y/nN  for  several  different  values  of  L  in  Fig.  7. 


10  Sampling  Procedure 

Finally,  we  consider  some  practical  considerations  concerning  the  computer  implementation  of  the  ILGA 
model.  Since  each  site  has  nL  bits,  and  therefore  2”^  possible  states,  and  since  the  most  popular  lattices 
with  the  requisite  isotropy  properties  have  n  =  6  and  n  =  24,  it  is  clear  that  the  brute-force  approach  in 
which  a  lookup  table  is  used  to  store  the  collision  outcome  states  will  not  be  feasible  for  L  much  greater 
than  unity. 

For  this  reason,  we  propose  another  sampling  scheme  for  the  outgoing  states.  Though  the  method  we 
propose  is  completely  general,  we  illustrate  it  for  the  two  dimensional  integer  lattice  gas  on  a  triangular  grid 
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Fractional  Noise 


(n  =  6).  Let  n  be  an  integer- valued  column  n- vector  whose  components  are  the  particle  occupation  numbers 
in  each  of  the  six  directions. 

Let  us  suppose  that  the  mass  and  the  two  components  of  momentum  are  the  only  conserved  quantities. 
Since  these  conserved  quantitites  are  linear  in  the  particle  occupation  numbers,  each  of  them  correspond  to 
a  row  vector,  whose  inner  product  with  n  yields  the  conserved  quantity  in  question.  Thus,  corresponding  to 
the  mass  we  have  the  row  vector 

qi  =  (  1  1  1  1  1  1  ) , 

corresponding  to  the  a;-momentum  (multiplied  by  a  factor  of  2),  we  have 

q2  =  (  2  1  -1  -2  -1  1  ) , 

and  corresponding  to  the  j/- momentum  (multiplied  by  a  factor  of  2/\/3),  we  have 

qa  =  (  0  1  1  0  -1  -1  )  . 

In  fact,  these  row  vectors  are  precisely  the  hydrodynamic  eigenvectors,  mentioned  in  our  derivation  of  the 
viscosity;  that  is, 

A(qi)"  =  A(q2)*  =  A(q3)*  =  o. 

It  is  clear  that  these  can  always  be  chosen  to  be  mutually  orthogonal,  without  loss  of  generality.  Using  the 
Gram-Schmidt  procedure,  it  is  then  possible  to  find  three  vectors  spanning  the  kinetic  subspace,  orthogonal 
to  the  above;  e.g., 

q4  =  (  2  -1  -1  2  -1  -1  ), 
q5  =  (  1  -1  1  -1  1  -1  ), 
and 

qe  =  (  0  1  -1  0  1  -1  ) . 

Now  the  collision  process  takes  state  n  to  state  n'.  Since  it  cannot  change  the  values  of  the  conserved 
quantities,  it  follows  that  the  difference  between  n'  and  n  must  be  a  linear  combination  of  kinetic  eigenvectors. 
That  is, 

n'  =  n  -b  a4qj  +  asqj  +  Oeq^, 

where  the  a’s  are  integer  constants,  and  where  the  superscript  T  denotes  “transpose.”  Thus,  writing  out 
components,  we  have 

/  ni  -f-  2a4  -h  Os  ) 
n2  —  0!4  —  0:5  +  t»6 
,  _  na  —  04  +  —  ae 

TI4  2a4  —  0/5 
ns  —  04  +  05  +  ae 
\  n6  —  a4  —  as  —  as 

Since  the  components  of  n'  must  all  be  between  0  and  2^  —  1,  inclusive,  we  derive  the  following  six  inequality 
constraints: 


0  <  ni  +  2a4  -b  as  <  2^  -  1 

0  <  n2  —  04  —  05  -b  ae  <  2^^  —  1 

0  <  na  —  04  -b  as  —  ae  <  2^  —  1 

0  ^  ^4  -b  2a4  —  as  ^  2^  —  1 

0  <  ns  —  a4  -b  as  -b  ae  <  2^  —  1 

0  <  ne  —  a4  —  as  —  oe  <  2^  —  1. 
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These  inequality  constraints  define  a  polytope  in  the  three  dimensional  space  of  allowed  values  of  0:4,  as  and 
as-  We  know  that  this  pol)d;ope  must  exist  and  contain  the  origin  in  that  space,  since  =  05  =  ae  —  0, 
corresponding  to  the  “trivial  collision”  in  which  the  occupation  numbers  do  not  change  their  values,  will 
always  satisfy  the  constraints. 

The  collision  process  is  then  specified  by  a  strategy  for  sampling  points  from  this  polytope.  One  viable 
strategy  which  certainly  satisfies  semidetailed  balance  is  to  sample  the  points  within  this  polytope  uniformly. 
It  is  possible,  though  tedious,  to  derive  a  closed-form  algorithm  to  do  this,  based  on  the  above  constraints. 
Alternatively,  with  some  loss  of  efficiency,  one  can  simply  bound  the  polytope  and  use  a  rejection  sampling 
scheme.  Details  of  this  procedure  will  be  provided  in  a  forthcoming  publication  [14]. 

11  Conclusions 

We  have  generalized  the  hydrodynamic  lattice  gas  model  to  include  integer  numbers  of  particles  moving  in 
each  direction  at  each  site.  We  have  presented  the  thermodynamics  and  kinetic  theory  of  this  generalized 
Integer  Lattice  Gas  (ILGA)  model,  including  closed-form  (or  parametric  algebraic)  equations  for  the  equi¬ 
librium  distribution  function,  the  entropy,  the  equation  of  state,  the  non-galilean  factor  in  the  inertial  term 
of  the  fluid  equations,  and  the  statistical  noise.  We  have  thereby  shown  that  the  ILGA  model  allows  for  the 
attainment  of  galilean  invariance,  and  a  reduction  in  the  kinematic  viscosity  and  the  statistical  noise.  In 
future  publications,  we  shall  show  that  this  generalization  also  allows  for  more  straightforward  inclusion  of 
interparticle  interactions  than  the  usual  binary  model. 
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A  The  Infinite  Integer  Limit 

To  consider  the  limit  of  infinite  integers,  T  -4  oo,  we  first  note  that  the  fractional  occupation  number 

Fi,{z)  1  / 


h{z)  = 


has  the  limiting  behavior 


2^-1  2^  -  1  \  1 


lim  fL{z)  =  0 

z~^Q 


1 


(A.23) 


lim  h(z)  =  i 

Z-fl 


1.  Referring  to  Figure  1,  we  note 
►  oo.  To  verify  this,  we  note  that 


lim  /l(z)  =  1 

^-♦OO 

for  all  L\  here  we  have  used  L’Hopital’s  rule  to  establish  the  result  for  z  - 
that  the  function  /i(^)  becomes  increasingly  like  a  step  at  z  =  1  as  T 
the  width  of  the  gradient  there  can  be  estimated  by 

fU^)  6 

2^  +  1’ 

which  clearly  goes  to  zero  as  L  -4  oo;  once  again  we  have  used  L’Hopital’s  rule  to  establish  this  result. 

The  approach  to  a  step  function  means  that  the  entire  range  of  fractional  occupation  numbers  is 
parametrized  by  values  of  z  within  order  2~^  from  1,  as  L  ->  oo.  That  being  the  case,  we  write 

y 


Z  —  Id" 


2^’ 


(A.24) 


where  y  is  a  new  parameter  of  order  unity.  Note  that  the  fractional  occupation  number  is  exactly  1/2  when 
y  =  0.  Inserting  Eq.  (A.24)  into  Eq.  (A.23),  we  can  now  take  the  limit  as  L  -t  oo  to  get 


lim  +  = 

L‘-+oo  \  2^/  1—6”*'  2/ 


(A.25) 


Next,  inserting  Eq.  (A.24)  into  Eqs.  (4.9)  and  (4.10),  and  taking  the  limit  as  L  -t  oo,  we  find  the  entropy 
excess. 


^lim  (l  +  ^)  =  In  (1  -  e-^')  -  -Iny+1. 


(A.26) 


Eqs.  (A.25)  and  (A.26)  constitute  parametric  algebraic  equations,  with  parameter  y,  yielding  AS^  as  a 
function  of  the  fractional  occupation  number  /l  as  L  — >  oo.  These  equations  were  used  to  produce  the 
shaded  curve  in  Fig.  2. 

Likewise,  inserting  Eq.  (A.24)  into  Eqs.  (9.22),  and  taking  the  limit  as  L  — >  oo,  we  find  the  fractional 
noise, 


lim  AFi  fl  +  =  J-  2~^  ■ 

L-*oo  \  2^/  \Jy‘  —  2ysi 


j/2  —  2  cosh  y  +  2 


(A. 27) 

2j/sinh3/ +  2coshy  —  2  ’  ' 

Eqs.  (A.25)  and  (A. 27)  constitute  parametric  algebraic  equations,  with  parameter  y,  yielding  ATi  as  a 
function  of  the  fractional  occupation  number  fi,  as  L  —¥  oo.  These  equations  were  used  to  produce  the 
shaded  curve  in  Fig.  7. 

Finally,  inserting  Eq.  (A.24)  into  Eqs.  (7.19),  and  taking  the  limit  as  L  — >  oo,  we  find  the  Gl  factor  for 
a  Bravais  lattice, 

1-  r  ,  y\  2+{y^  +  2y-8)ey+{y^-6y+12)e^y+{y^-y^  +  6y-8)e^y-2(y-l)e^y 

L-^oo  2^)  (1  -  (y2  2)ey  + 

(A.28) 

Eqs.  (A.25)  and  (A.28)  constitute  parametric  algebraic  equations,  with  parameter  y,  yielding  Gl  as  a  function 
of  the  fractional  occupation  number  fi  as  L  oo.  These  equations  were  used  to  produce  the  shaded  curve 
in  Fig.  3. 
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